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ABSTRACT 

We test a method of estimating the power spectrum of turbulence in molecular clouds based on the 
comparison of power spectra of integrated intensity maps and single- velocity-channel maps, suggested 
by Lazarian and Pogosyan. We use synthetic ^'^CO data from non-LTE radiative transfer calculations 
based on density and velocity fields of a simulation of supersonic hydrodynamic turbulence. We 
find that the method yields the correct power spectrum with good accuracy. We then apply the 
method to the Five College Radio Astronomy Observatory ^^CO map of the Perseus region, from the 
COMPLETE website. We find a power law power spectrum with slope (3 = 1.81 ± 0.10. The values 
of /3 as a function of velocity resolution are also confirmed using the lower resolution map of the same 
region obtained with the AT&T Bell Laboratories antenna. Because of its small uncertainty, this 
result provides a useful constraint for numerical codes used to simulate molecular cloud turbulence. 
Subject headings: ISM: clouds - ISM: kinematics and dynamics ~ ISM: structure 



1. INTRODUCTION 

The large scale dynamics of interstellar clouds is char- 
acterized by random supersonic motions with a very large 
Reynolds number, meaning that inertia is much larger 
than viscous forces. This interstellar turbulence is a 
dominant transport mechanism in many astrophysical 
processes and plays an important role in the fragmen- 
tation of star- forming clouds. Random supersonic flows 
in molecular clouds result in a complex network of highly 
radiative shocks causing very large density contrasts and 
shaping the clouds into the observed self-similar filamen- 
tary structure (Nordlund & Padoan 2003). 

Our understanding of turbulence is limited by the lack 
of general analytical solutions of the Navier-Stokes equa- 
tion. Statistical properties of turbulence, essential to 
modeling many astrophysical processes, are therefore de- 
rived almost entirely from numerical simulations. In the 
case of incompressible hydrodynamic turbulence, numer- 
ical simulations can be compared with laboratory exper- 
iments. This is not possible for the isothermal, super- 
sonic, magneto-hydrodynamic (MHD) turbulence regime 
of molecular clouds. The best approach to validate nu- 
merical models is in this case a statistical comparison 
of observational data with synthetic data from numerical 
simulations. 

Padoan et al. (2004) found that the exponents of the 
velocity structure functions of compressible and super- 
Alfvenic turbulence follow a generalized She-Leveque 
scaling (Boldyrev 2002), depending only on the rms Mach 
number of the flow. Based on that scaling formula, if 
one of the exponents is known (for example the second 
order that corresponds to the velocity power spectrum) 
all the others can be derived, up to some order. While 
the scaling formula is well constrained numerically, its 
normalization, given by the actual value of one of the 
exponents, is harder to measure and may depend on the 
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numerical method used to simulate the turbulence. 

Numerical simulations of supersonic hydrodynamic 
(HD) and magneto-hydrodynamic (MHD) turbulence 
have yielded a range of values of the slope of the velocity 
power spectrum, from the Kolmogorov value of /3 = 5/3, 
to the Burgers value of /3 = 2 (Frisch & Bee 2001), and 
beyond. The main problems with the numerical esti- 
mate of the power spectrum are i) The limited extent 
in wavcnumbers of the inertial range of turbulence, or 
even its complete absence in the case of low resolution 
(or highly dissipative) simulations; ii) The emergence of 
the bottleneck effect (e.g. Falkovich 1994; Dobler et al. 
2003; Haugen & Brandenburg 2004) as soon as the nu- 
merical resolution is large enough to generate an inertial 
range; iii) The dependence of the power spectrum on 
the numerical schemes used to stabilize the shocks; iv) 
The dependence of the numerical resolution necessary for 
convergence on the numerical method. Given these dif- 
ficulties, reliable measurements from interstellar clouds 
provide useful constraints to validate numerical models 
and to improve our knowledge of supersonic turbulence. 

Lazarian & Pogosyan (2000) have demonstrated that 
the exponent of the velocity power spectrum can in prin- 
ciple be derived from spectral maps of emission lines by 
comparing the power spectrum of integrated intensity 
with the power spectrum of single- velocity-channel inten- 
sity. Their method was tested with numerical simulations 
of turbulent flows, where the velocity and density flelds 
had to be modified to generate power law power spec- 
tra, due to the limited numerical resolution (Esquivel 
et al. 2003). In this Letter we present a new test of 
their method. Thanks to the high numerical resolution of 
our simulation (1, 024'^ computational zones), our density 
and velocity power spectra exhibit extended power laws, 
and no modifications to the original turbulent fields are 
required. Furthermore, the method is tested by comput- 
ing synthetic CO emission lines with a non-LTE radia- 
tive transfer code. We find the method allows to retrieve 
the power law exponent of the velocity power spectrum 
with good accuracy, using the J=l-0 ^^CO line. We then 
apply the method to the Five College Radioastronomy 
Observatory (FCRAO) survey of the Perseus molecular 
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Simulation: j8=1.83±0.17 




Fig. 1. — Compensated power spectrum of integrated intensity 
(bottom curve) and average of the compensated power spectra of 
single- velocity-channel intensity maps (top curve) from the numer- 
ical simulation. The dotted plot (middle curve) is the 3D power 
spectrum of velocity and the dashed line shows the /3 = 1.83 slope 
predicted by the method. 



cloud complex (Ridge at al. 2006), publicly available 
from the COMPLETE website, and estimate the power 
spectrum exponent /? = 1.81 ± 0.10. The result is con- 
firmed by the analysis of the lower resolution map of the 
same region obtained with the AT&T Bell Laboratories 
antenna (see Padoan et al. 1999). 

2. SYNTHETIC SPECTRAL MAPS 

Synthetic maps of the J=l-0 line of ^^CO are com- 
puted with a new non-LTE radiative transfer code that 
has been extensively tested against the older Monte Carlo 
code (Juvela 1997). As a model of the density and ve- 
locity fields in interstellar clouds we use the results of a 
supersonic simulation of HD turbulence with rms Mach 
number M = 6. Our purpose is not to simulate the con- 
ditions found in a specific molecular cloud, but rather to 
test a general method of estimating the velocity power 
spectrum. We will study the effect of different parameter 
values elsewhere. Because the method is based on an an- 
alytical derivation by Lazarian & Pogosyan (2000) that 
neglects correlations between density and velocity, we ex- 
pect it to work equally well for a vast range of turbulence 
parameters. 

The simulations are carried out with the Enzo code, 
developed at the Laboratory for Computational Astro- 
physics by Bryan, Norman and collaborators (Norman 
& Bryan 1999). Enzo is a public domain E ulerian 
grid-based code (see http://cosmos.ucsd.edu/enzo/) that 
adopts the Piecewise Parabolic Method (PPM) of Colella 
& Woodward (1984). We use an isothermal equation 
of state, periodic boundary conditions, initially uniform 
density and initial random large scale velocity. The tur- 
bulence is forced in Fourier space only in the wavenumber 
range 1 < k < 2, where k = 1 corresponds to the size of 
the computational domain that contains 1 , 024"^ compu- 
tational zones (for details see Kritsuk et al. 2006). 

The radiative transfer calculations assume a box size 
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Fig. 2. — Convergence plot for the exponent of the velocity power 
spectrum as a function of the channel width. 



of 5 pc, a mean density of 10^ cm~^, a mean kinetic 
temperature of 10 K, an rms Mach number M = 6 (con- 
sistent with the HD simulation) and a uniform abun- 
dance of ^"^CO molecules. These values were chosen as 
a generic reference model, not tailored to the Perseus 
molecular cloud complex. The original density and ve- 
locity data-cubes are re-sampled to a resolution of 256^ 
zones. This has three advantages: i) It leaves density and 
velocity fields with power spectra that are power laws 
almost up to the new Nyquist frequency; ii) It speeds 
up enormously the radiative transfer calculations; iii) 
It generates a map of synthetic spectra with a range of 
scales comparable to that of the rebinned FCRAO map 
of Perseus (see below) . 

The result of the radiative transfer calculations is a 
map of spectral profiles of the line intensity, T(w,x), 
over 280 velocity channels, v, and across all positions, 
X, within the map of 256^ positions. The sum of the 
line intensity of all velocity channels, multiplied by the 
channel width, Av = 0.025 km/s, gives the integrated in- 
tensity at each position, /(x) = Ei,T(u,x) Av. The map 
/(x) may be used to derive a rough estimate of the col- 
umn density. We do not address the issue of converting 
integrated intensities into accurate estimates of column 
density, because we derive the velocity power spectrum 
directly from the line intensity data. 

3. THE METHOD 

We call Pi(fc) the two-dimensional spatial power spec- 
trum of the integrated intensity map, /(x), and look for 
a power law fit such that, Pi{k) oc The integrated 

map is the sum of all the single-velocity-channel maps. 
We then call Pt{v, k) the two-dimensional spatial power 
spectrum of the individual channel maps, r(u,x), where 
we have left the velocity dependence to stress that there 
are as many of these maps and power spectra as there 
are velocity channels. We look for a power law fit to 
the average of these power spectra as well, such that 
Prj.{k) = {PT{v,k))^ (X fc-Z^T. Finally, we assume the 
three dimensional spatial power spectrum of the velocity 
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Perseus: |8=1.81±0.10 
/St= 1.59 ±0.01 
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Fig. 3. — Same as in Figure [T] but for the FCRAO survey of 
Perseus (Ridge et al. 2006). 



field is also a power law (a reasonable assumption for 
the cloud turbulence, given the huge Reynolds number), 

The result of Lazarian & Pogosyan (2000) can be ex- 
pressed as /3 = 1 + 2(/3i - /3t), if A < 2, and /3 = 
1 + 2(2 - /3t), if A > 2. In this Letter we test their 
result only for the case of /3i < 2. The single- velocity- 
channel power spectra are shallower than the integrated 
intensity power spectrum, /?t < A, because they contain 
small-scale structure originating from both the density 
and the velocity fields. It is this difference that makes it 
possible to derive the velocity power spectrum exponent, 
based on the above formula. 

In this Letter we refer to power spectra as the total 
power within shells of wavenumber between k and k + dk. 
In Lazarian & Pogosyan (2000), the power spectra are 
defined as the average power within shells of wavenumber 
between k and k + dk. Furthermore, we do not absorb 
the negative sign in the definition of /?/ as they do. Our 
exponents are related to their n and m exponents in the 
following way: /3/ = 1 — n, /3 = 1 -I- m. In our convention, 
the velocity power spectrum corresponds to the usual 
tmbulent energy spectrum, for example k~^/'^ for the 
Kolmogorov case. 

Fig. [T] shows Pi{k) (bottom plot) and Prik) (top plot) 
for the simulated data-cube "observed" along the x di- 
rection. The error bars are one standard deviation 
above and below the mean values and depend on the 
statistical sample size (number of Fourier modes inside 
each wavenumber shell), so they decrease with increasing 
wavenumber. The linear least square fits, plotted as solid 
lines, are computed in the wavenumber range I < k < 30. 
Based on the estimated exponents of these two power 
spectra, we obtain P = 1 + 2(/3i - /3t) = 1.83 ± 0.17. 
Values of (3 from synthetic maps of the data-cube ob- 
served in different directions fall within the estimated 1- 
(T uncertainty. The value of the velocity power spectrum 
exponent, computed directly from the original three di- 
mensional velocity field, is (3 = 1.8 ± 0.1 (dotted plot in 
Fig. [T]), showing that the method retrieves the correct 



exponent. 

In order to verify the dependence of (3 on the velocity 
channel width, we have applied the method to synthetic 
maps at different velocity resolutions. The result, plotted 
in Fig. O shows that the value of /3 is completely con- 
verged only for a channel width of the order of the ther- 
mal width. This was to be expected, because all velocity 
fluctuations above the thermal width may in principle 
affect the velocity power spectra of the single-velocity- 
channel maps. 

This method is based on an analytical derivation by 
Lazarian & Pogosyan (2000) that neglects the correla- 
tions of density and velocity in turbulent flows. We in- 
terpret this test as a confirmation of the validity of their 
method in the case of /3i < 2, rather than as an empirical 
calibration of the value of /3 based on /3i and Pt- The un- 
certainty of the method, when applied to observational 
data, is then determined by the error bars of the obser- 
vational power spectra, independently of the uncertainty 
of our numerical test. The final error bar is dominated 
by the uncertainty in the power spectrum exponent of 
the integrated intensity, because the uncertainty in the 
single- velocity-channel spectrum is reduced by averaging 
the power spectra of many velocity channels. 

4. THE POWER SPECTRUM OF PERSEUS 

We have appHed the method to the J=l-0 ^'^CO sur- 
vey of the Perseus molecular cloud complex carried out 
with the FCRAO 14 m antenna by Ridge et al. (2006). 
The grid spacing of the survey is 23" , and the beam size 
46". The velocity-channel size is 0.06 km/s. The power 
spectra we compute are corrected for the effect of beam 
and noise, by simply dividing by the power spectrum of 
a gaussian beam, and subtracting the power spectrum of 
the noise. Spatial correlations in the noise arising from 
the "On-the-Fly" mapping mode are neglected. At the 
largest wavenumbers, the power spectra are sensitive to 
the noise subtraction, and realistic error bars account- 
ing for that would make such wavenumbers essentially 
useless for estimating the power spectra. We therefore 
prefer to regrid the map to a resolution of 92" , which 
has the advantage of increasing the signal-to-noise by a 
factor of 4. 

The power spectra are shown in Fig. [3l The least 
square fits are computed in the range 5 < A: < 80, and 
yield /3 — 1.81 ± 0.10. Although this range is less ex- 
tended than that used with the synthetic data, it includes 
larger wavenumbers (between fc = 30 and k — 80) than 
the synthetic fit (the synthetic power spectra are still 
affected by numerical dissipation at large wavenumbers 
even after rebinning from 1,024'^ to 256"^). This reduces 
the final uncertainty, because the statistical sample size 
is much larger at larger wavenumbers (oc k dk in two di- 
mensions). As a result, the 1-cr uncertainty of f3 is 9% 
for the synthetic data, and 5% for the observations. 

The value of the slope of the projected density power 
spectrum, /3i — 1.99 ±0.05 is similar to values previously 
found in different regions observed in HI (e.g. Green 
1993; Stanimirovich et al. 1999) and CO (e.g. Bensch, 
Stutzki, Ossenkopf 2001; Padoan et al. 2004a). No- 
tice that in those previous works, with the exception of 
Padoan et al. (2004a), the power spectrum is not inte- 
grated over wavenumber shells, so its slope is equivalent 
to Pi + 1. Furthermore, the power spectrum of the ^^CO 
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integrated intensity is slightly steeper than the gas den- 
sity power spectrum due to radiative transfer effects. Ac- 
counting for such effects, the slope of the projected gas 
density power spectrum in the Perseus region was esti- 
mated to be consistent with that of super- Alfcnic turbu- 
lence simulations (Padoan et al. 2004a). 

The value of /3 estimated for the Perseus region as a 
function of the velocity resolution is shown in Fig. [2] (solid 
line). The channel width of 0.06 km/s is close to the 
thermal line width and the value of (3 seems to be al- 
most converged, at least within its 1-cr uncertainty. As 
an independent test, we have applied the method also 
to the AT&T Bell Laboratories map of the same region 
(see Padoan et al. 1999). The Bell Laboratories 7 m an- 
tenna has a beam twice the size of the FCRAO 14 m an- 
tenna. We did not regrid this map to a lower resolution, 
so the map resolution and size are comparable to those 
of the FCRAO map, but its velocity resolution is only 
0.273 km/s, so it is not expected to yield a converged 
value of (3. The values of /3 for the Bell Laboratories 
map is shown in Fig. [5] as a dotted line, showing very 
good agreement with the FCRAO result, well within the 
estimated 1-cr uncertainty. 

5. DISCUSSION AND CONCLUSIONS 

This result has very interesting implications. First, 
numerical simulations giving power spectra with slope 
significantly larger than /3 = 1.8 ± 0.1 may be ruled out 
as correct descriptions of molecular cloud turbulence (at 
least for the Perseus region). Burgers exponent, /3 = 2.0, 
is 2a larger than the Perseus exponent (assuming this 
is converged as a function of velocity resolution). The 
slope of the power spectra of the SPH simulations in 
Ballesteros-Paredes et al. (2006), (i k, 2.7 in the case of 
a turbulence rms Mach number « 6, is 9(T larger than the 
present estimate, and their grid based simulations have 
[3 « 2.2, 4cr too large. Second, we can now derive the 
absolute values of the velocity structure function expo- 
nents in Perseus. Padoan et al. (2004b) have determined 
numerically the relative values of those exponents, so the 



knowledge of one of them, for example the second order 
exponent given by — 1, allows to determine the absolute 
values of all the others. 

A different method of estimating the scaling of the tur- 
bulence from molecular clouds surveys was developed by 
Brunt & Heyer (2002a, b), based on the principle compo- 
nent analysis (PCA). They analyzed 23 molecular clouds 
in the outer Galaxy and estimated a value fi = 2.2 ± 0.3, 
if the structure function exponents of order p are as- 
sumed to scale linearly as p/3. This value is significantly 
larger than the one we estimate in Perseus. However, the 
PCA method is dependent on a calibration with numer- 
ical simulations. Based on such a calibration, it appears 
that the PCA method estimates the exponent of the ve- 
locity structure fimctions of order p — 0.5 or lower (Brunt 
et al. 2003). Taking this into consideration, the result of 
Brunt & Heyer (2002b) would be roughly consistent with 
ours, if a very intermittent scaling of the structure func- 
tion is assumed, consistent with numerical simulations of 
supersonic turbulence. 

We have shown in this Letter that the method for es- 
timating the velocity power spectrum slope proposed by 
Lazarian and Pogosyan (2000) works well in the case of 
/3i < 2, and for a velocity resolution not much larger 
than the thermal line width. However, regions with lower 
turbulence Mach number than Perseus, for example the 
Taurus region, may have steeper density power spectra, 
and hence (3i> 2. Such regions will be studied in a future 
work, where the method will be tested also with numer- 
ical simulations generating synthetic data with /3i > 2. 
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